#########################################################
# Name: Trevor Osaki                                    #
# Date: Oct 11, 2023                                    #
# Purpose: Test some hypotheses (Appendix 7.3)          #               
#########################################################

# Here, we should be focusing on Group 2 individuals ("utilitarians")
# We want to make an inference about alpha of paper
# In short, we do this by working with the Wald Statistic
# Other file alpha_appendix_A7_4.R is arranged similarly

### Set up some parameters ###
# These are taken from results (from "4_tables.do")
# Mean fairness evaluation toward Anti-Black Discrimination (Stage 2)
b2 = -0.8869565;
# Mean fairness evaluation toward Anti-White Discrimination (Stage 2)
w2 = -0.7083333;
# Mean fairness evaluation toward Anti-Black Discrimination (Stage 1)
b1 = -1.407643;
# Mean fairness evaluation toward Anti-White Discrimination (Stage 1)
w1 = -0.377193;

# Set up vector of means
est = c(b2,w2,b1,w1);

# Define the variance and covariances
v_b2 = 0.02063231;
v_w2 = 0.02604875;
v_b1 = 0.01208239;
v_w1 = 0.01373570;

v_b2_w2 = 0;
v_b2_b1 = 0;
v_b2_w1 = 0.00772206;
v_w2_b1 = 0.00996603;
v_w2_w1 = 0;
v_b1_w1 = 0;

# Covariance matrix setup
var = c(v_b2, v_b2_w2, v_b2_b1, v_b2_w1, v_b2_w2, v_w2, v_w2_b1, v_w2_w1, v_b2_b1, v_w2_b1, v_b1,v_b1_w1,v_b2_w1, v_w2_w1,v_b1_w1,v_w1);

# Sequence/vector of possible alpha values
alpha = seq(-1,2,0.001);

# Number of restrictions
q = 1;

# Number of parameters
k = 4;

### Set up matrices ### 
# All items will be stored in "test"
test  = data.frame(matrix(0, nrow = length(alpha),ncol = 5))
colnames(test) <- c("alpha", "wald_WB", "p_WB","wald_BW","p_BW")
# Vector of means
theta = matrix(est,nrow = k,ncol = 1,byrow = TRUE);
# Vector of scalars to use in calculating the Wald statistic
r     = matrix(0,nrow = q,ncol = 1);
# Covariance matrix
V     = matrix(var,nrow = k,ncol = k,byrow = TRUE);

for(i in 1:length(alpha)){
  
  test$alpha[i] = alpha[i];
  
  # For calculating Wald statistic among W-B switchers
  scalars1 = c(1,0,-alpha[i],-(1-alpha[i]));
  # For calculating Wald statistic among B-W switchers
  scalars2 = c(0,1,-(1-alpha[i]),-alpha[i]);
  
  R1 = matrix(scalars1,nrow = q,ncol = k);
  R2 = matrix(scalars2,nrow = q,ncol = k);
  
  R1t = t(R1);
  R2t = t(R2);
  
  B1 = solve(R1%*%V%*%R1t);
  B2 = solve(R2%*%V%*%R2t);
  
  A1 = R1%*%theta - r;
  A2 = R2%*%theta - r;
  
  A1t = t(A1);
  A2t = t(A2);
  
  # Calculate Wald Statistic for W-B switchers
  wald1 = (A1t%*%B1%*%A1);
  # ...and B-W switchers
  wald2 = (A2t%*%B2%*%A2);
  
  test$wald_WB[i] = wald1;
  p_wb = 1-pchisq(wald1,df = q);
  test$p_WB[i] = p_wb;

  test$wald_BW[i] = wald2;
  p_bw = 1-pchisq(wald2,df = q);
  test$p_BW[i] = p_bw;
}

# Only keep "test"
rm(list=setdiff(ls(), "test"))

  # In "test,"
  # Col 1 contains possible values of alpha
  # Col 2 contains the Wald statistic for W-B switchers
  # Col 3 contains p-values
  # Col 4 contains the Wald statistic for B-W switchers
  # Col 5 contains p-values


# 1. Find the alpha that yields the highest p-value
cat("alpha=", test[which.max(test$p_WB),]$alpha, "for the White-to-Black switchers.") # 0.495
cat("alpha=", test[which.max(test$p_BW),]$alpha, "for the Black-to-White switchers.") # 0.679

# 2. Find p-values when alpha = 0, 0.5, 1
selectedAlpha <- test[test$alpha == 0.000 | test$alpha == 1.000 | test$alpha == 0.500, 
                 c('alpha', 'p_WB', 'p_BW')]
selectedAlpha$p_WB <- round(selectedAlpha$p_WB, 3)
selectedAlpha$p_BW <- round(selectedAlpha$p_BW, 3)
selectedAlpha

# 3. Back out CI for W-B switchers via calculating p-value
  # This uses the fact that the Wald is distributed chi2(q).
  # The alphas with p_wb~0.05 form the approx. bounds
p_WB_indice_low <- which(diff(sign(test$p_WB - 0.05)) == 2)
p_WB_indice_up <- which(diff(sign(test$p_WB - 0.05)) == -2)
p_WB_indices <- c(p_WB_indice_low, p_WB_indice_low+1, p_WB_indice_up, p_WB_indice_up+1)
test[p_WB_indices, c('alpha', 'p_WB')]
  # The CI should be [0.242,0.801].


# 4. Back out CI for B-W switchers via calculating p-value
  # The alphas with p_bw~0.05 form the approx. bounds
p_BW_indice_low <- which(diff(sign(test$p_BW - 0.05)) == 2)
p_BW_indice_up <- which(diff(sign(test$p_BW - 0.05)) == -2)
p_BW_indices <- c(p_BW_indice_low, p_BW_indice_low+1, p_BW_indice_up, p_BW_indice_up+1)
test[p_BW_indices, c('alpha', 'p_BW')]
  # The CI should be [0.404,1.076]

